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Abstract 

''sT "By implementing a recent technique for the determination of stochastic eigendirections of two coupled stochastic variables, we 
.investigate the evolution of fluctuations of NO2 concentrations at two monitoring stations in the city of Lisbon, Portugal. We analyze 
,— i"the stochastic part of the measurements recorded at the monitoring stations by means of a method where the two concentrations 
Ph .are considered as stochastic variables evolving according to a system of coupled stochastic differential equations. Analysis of their 
^ "structure allows for transforming the set of measured variables to a set of derived variables, one of them with reduced stochasticity. 
d .For the specific case of NO2 concentration measures, the set of derived variables are well approximated by a global rotation of the 
"original set of measured variables. We conclude that the stochastic sources at each station are independent from each other and 
'"O .typically have amplitudes of the order of the deterministic contributions. Such findings show significant limitations when predicting 
C/j 'such quantities. Still, we briefly discuss how predictive power can be increased in general in the light of our methods. 
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Figure 1: NO2 measurement stations in the region of Lisbon (Portugal) at tlie 
Soutliwestern coast of Europe. In this paper we focus on the set of measure- 
ments taken at the stations of Chelas and Avenida da Liberdade with approxi- 
mately 10^ data points each extracted in the period between 1995 and 2006. 



1. Introduction 

The industrial and urban development during the last decades 
has led to a general decrease of air quality, drastically affecting 
urban environmental and human life quaUty. Although, accord- 
ing to the European Environment Agency reportfH, air quality 
has improved in general during the last years, this enhancement 
was not significant enough to ensure good air quality in all ur- 
ban areas. One of the pollutants with negative impact on health 
and environment is NO2. Anthropogenic NO2 is mainly emitted 
by vehicles and industrial processes. NO2 has not only severe 
effects on health causing e.g. respiratory and cardiovascular dis- 
eases, it also affects the environmentyl as nitrogen deposition 
leads to eutrophicationQ. A better understanding of the mech- 
anisms that influence production, transport, and decomposition 
of NO2 is therefore important. Previous studies revealed that 
temperature, wind speed and direction, relative humidity, cloud 
cover, dew point temperature, sea level pressure, precipitation. 



and mixing layer height are relevant meteorological variables 
to model the concentrations of air pollutantsilllllHIl. In 
particular, approaches that deal with the evolution of the NO2 
concentration at individual city locations are important for fore- 
casting the air quality of urban regions. 

Recently a frame work ||^[l3l for analyzing measurements on 
complex systems was introduced, aiming for a quantitative es- 
timation of drift and diffusion functions from measured data. 
These functions can be identified with the deterministic and 
stochastic contributions to the dynamics, respectively, and give 
a considerable insight into the underlying systems. The frame- 
work was already successfully applied for instance to describe 
turbulent flows||9|] and the evolution of climate indices fll lll29ll . 
stock market indices 1 121. and oil prices joj- At the same time. 



the basic method has been refined in particular with respect to 
the impact of finite sampling effects lll4ll5ll . the impact of mea- 
surement noise 1 3 1/7, 18], and the role of local eigendirections 
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of the diffusion matrices lfl9ll . 

In this paper, we aim to apply some recent methods for de- 
riving variables with reduced stochastic fluctuations to empiri- 
cal data. Namely, we adapt this framework for analyzing mea- 
surements of NO2 concentrations in the metropolitan region of 
Lisbon, Portugal (see Fig. [1]), taken over several years. We ar- 
gue that the temporal fluctuations of these concentrations re- 
sult from two independent contributions: one periodic and one 
stochastic. The periodic part describes daily, weekly, seasonal 
and yearly variations of the concentration, which is an accepted 
and well-studied resultj^l. The stochastic contribution can be 
modelled through a stochastic differential equation ll23ll having 
two terms, one drift forcing (deterministic) and one diffusive 
fluctuation (stochastic). 

When addressing stochastic higher-dimensional systems it is 
typically difficult to identify the variables most relevant for a 
proper description of the system's evolution. In geophysical 
applications, the reduction of the full set of variables to only a 
few variables often is achieved bymeans of the so-called Princi- 
pal Component Analysis (PCA') i2lll or other standard reduction 
methods, such as stepwise regression or ARIMAlS] . However, 
the inherent fluctuations are not so commonly investigated. 

In this paper we will apply a recent method for reconstruct- 
ing the phase space of two stochastic variables, which evolve 
according to a set of two coupled stochastic equations defined 
through drift vectors and diffusion matrices |19ll . The method is 
based on the eigenvalues of the diffusion matrices, from which 
it is possible to derive a path in phase space through which the 
deterministic contribution is enhanced. This technique implies 
a transformation of variables and allows for the investigation 
of the minimal number of independent sources of stochastic 
forcing in the system. In particular, a rather small eigenvalue 
of the diffusion matrix, compared to the average value of all 
the other, corresponds to one eigendirection in which stochastic 
fluctuations may be neglected, reducing the number of stochas- 
tic variables taken for describing the system's evolution. On the 
contrary, having all eigenvalues of the same order of magnitude 
means that the number of independent stochastic forces equals 
the number of variables. Moreover, as we will see, a direct in- 
spection of the diffusion functions enables one to ascertain if 
the stochastic contributions, one for each variable, are coupled 
among them or not. Therefore, we argue that the diagonaliza- 
tion of the diffusion matrices gives insight into the system. 

We start in Sec.|2]by describing the properties and the prepa- 
ration of the data set. Consecutively, in Sec.[3]and|4l the mod- 
eling of the time series as a Langevin process is carried out and 
its transformation to a new coordinate system are described in 
Sees. |5] and |6] respectively. In Sec. Qwe discuss the perfor- 
mance of the transformation of the coordinates obtained by our 
approach compared to other techniques commonly used for sta- 
tistical analysis of measured data. Section [8] closes this Letter 
with a general summary and ideas on the interpretation of the 
transformed time series with respect to the underlying environ- 
mental processes. 
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Figure 2: (a) Time series of tlie NO2 concentration at tlie station of Chelas, 
before detrending according to Eq. (2), and (b) a zoom-in of these "raw" time 
series yi compared to tlie detrended series .vi , wliich takes averages of 52-weeks 
periods, and then a second detrend with daily averages. Vertical offset of same 
plots are done for clarity. For the station at Avenida da Liberdade similar fea- 
tures are found (not shown). 

2. NO2 measurements in Lisbon 

In this section we briefly describe the sets of data analyzed in 
this paper as well as its preparation for analyzing the stochastic 
components of the measurements. 

The data set covers hourly measurements of NO2 concentra- 
tion, taken at 22 stations in the urban center of Lisbon recorded 
from of 1995 to 2006. For this study we choose the data 
from 1995 to 2005 for the monitoring stations at Chelas and 
at Avenida da Liberdade. These stations are located at a dis- 
tance of ~ 4 km from each other, see Fig.[T] In the following, 
the NO2 concentrations at the stations of Chelas and Avenida 
da Liberdade will be designated as yi (f) and y2{t), respectively, 
omitting the temporal dependency when not necessary. 

Increments in time are always of 1 hour. Each of the data 
sets contains 10^ measurement points approximately, including 
some periods of incomplete or erroneous measurements that are 
disregarded for our analysis. In the case of the chosen stations, 
the series of measurements yi and y2 contain 3726 and 4548 
instances of measurement errors, respectively. 

The concentration of NO2 is strongly driven by daily, weekly, 
monthly and yearly anthropogenic routines, and also by peri- 



3 



odic atmospheric processes. For instance the rush hours on 
working days have an almost immediate impact on the NO2 
concentration, and thus, on air quahty. The 24 hours and one 
week cycles are both traffic related and mirror daily and weekly 
cycles. The measurements of NO2 are therefore influenced by 
different periodic forcings and, since we are interested in the 
fluctuations of NO2 concentrations, the periodic behavior must 
be first detrended. The detrended series for yi and y2, repre- 
sented below as xi and X2 respectively, are obtained as follows. 

One first partitions the data in segments of length A^, which 
we suppose to be a multiple of relevant periodic fluctuations in 
the data set. As a second step, a mean segment is calculated by 
averaging measurements with the same position in the segment 
over the entire data set according to 



Uiin) s {yi{t)\t = n + niN, m = 0, 1, . . .) 



(1) 



for n = 0, 1, - 1. The detrended data set x,- is then calcu- 
lated by subtracting the respective values of the mean segment 
from the measured data, 



Xi{t) = y,(f) - Ni(t mod N) . 



(2) 



for t - I, . . .,T with T > N. If T is the size of the data set, 
our simulations have shown that averages over N = 52 weeks 
is the best choice for the entire data set, to take into account 
all known periodicities mentioned above. With this de trending 
method, some periodicities with variable phase remain. To filter 
also these periodicities, a second detrending with N = I day is 
then performed on consecutive periods of T = 14 days. 

Figure|2^ shows the original data yi for the station of Chelas. 
A zoom-in of a small time interval is plot in Fig. |2j3 together 
with the corresponding detrended data xi . From now on, if not 
stated explicitly otherwise, we will only consider the detrended 
time series xi and X2- Next describe their characteristics by 
means of a stochastic process. 

3. Modeling stochasticity in series of NO2 concentrations: 
Langevin processes 

The detrended series x, in Eq. (|2]i reflect the remaining 
stochastic components of the measurements at the respective 
stations of Chelas and Avenida da Liberdade. In this section 
we assume that, with two variables, the stochastic process is 
modeled by a system of two coupled Langevin equations, con- 
taining a deterministic and a stochastic part, described through 
a drift vector and a diffusion matrix, respectively. 

For the general case of a /iT-dimensional state vector X = 
(xi, xk), the Ito-Langevin equations describing the evolution 
of a particular trajectory in time read 1,2 3^ i241 : 



^ -h(X) + g(X)r(f), 

at 



(3) 



where F = (Fi , . . . , Tk) is a set of K independent stochastic 
forces with Gaussian distribution fulfilling 
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Figure 3: First conditional moments M''' for (a) the original series and (b) the 
detrended series, with different NO2 concentrations (xi.xi) at each one of the 
two stations (see legend). The corresponding second conditional moments M'^* 
are shown in (c) and (d), respectively. These moments are computed according 
to Eqs. and )9bt . While the original data presents oscillations beyond a 
given time interval ~ 5. the detrended time series does not (see text). The 
value of the corresponding Kramers-Moyal coefficient at the value of (xi , X2) 
chosen is given by Eq. (8) for the lowest value of r, i.e. one. 



The two terms on the right hand side of Eq. ^ include both the 
deterministic contribution, h = {/i,}, and the stochastic contri- 
bution, g = {gij}- The deterministic contribution describes the 
physical forces which drive the system, while functions g ac- 
count for the ampUtudes of the different sources of fluctuations 

The coefficients h and g are directly related to the drift vec- 
tors and diffusion matrices 1231] 



D\'\X) = /i,(X) 

K 

Df^iX) = ^g,-,(X)g,-,(X) 



(5) 
(6) 



k=\ 



for i,i - \,. . .,K, describing the evolution of the joint proba- 
bility density function (PDF) /(X, t) by means of the Fokker- 
Planck equation ||23 , 24 1 : 



K K 

■EE 

k=\ m=\ 



dxidxir, 



D['(X)f{X,t). 



(7) 



As done previously in other contexts^ [U E [B Q E 
1711 . the drift vector and the diffusion matrix can be derived di- 
rectly from the data. 

Statistically, the drift and diffusion coefficients coefficients 
Z)*'^ and Z)^^' are defined as 



D«(X).limiM!!(^, 

T^or kl 



(8) 




Figure 4: For the detrended series we plot (a-b) both components of the drift vector h = (/ii, /12) and (c-e) the components of diffusion matrix D*^' = {-D*^'}. The 

2^21 ■ 



coiTesponding fitted surfaces (black) are vertically offset for claiity. Since D*^' is symmetric (see text) one has Z)p* : 



with the first and second conditional moments given by 



M<'\x,T) = <y,(f + T)-y,(OIY(f) = x> 

■(Yj(t + T)-Yj(t)mt)^x)., 



M<f(X,T) 



(9a) 



(9b) 



Here (•|Y(f) = X) symbolizes conditional averaging over all 
events that fulfill the condition Y(f) = X. 

To determine the underlying Langevin equations, defined in 
Eq. O), one additionally needs to solve Eqs. (|5]l and In par- 
ticular, the calculation of matrices g from the diffusion matrices 
requires to solve D*'^' = gg^, e.g. by means of diagonalization, 
^diag ~ PD^^'P"', with P the orthogonal matrix of eigenvectors 

of D<^>. The family of solutions is given by g = P^ yjo^PO, 

where O is an arbitrary orthogonal matrix, obeying OO^ = 1 . 
The matrices D*^' are symmetric and positive semi-definite with 
all their eigenvalues real and non-negative (see Eq. (|9b]i). and 

therefore -^D^j^ is well-defined. For any choice of O the anal- 
ysis below does not change, and therefore we choose for sim- 
plicity O as the identity matrix. 

The computation of the conditional moments is based on 
their statistical r-dependence for small ITtIi . Previous 
works showed that Eqs. ( |9al ) and ( |9b] i are an operational defini- 
tion of the conditional moments that can easily be implemented 
for the direct estimation of the drift and diffusion coefficients 
from the data lfioi [itIi . In some practical situations, the limit 
in Eq. (O can be approximated by the slope of a linear fit of 
the corresponding conditional moments at small t. When such 
linear fit is not possible, an alternative estimate is to consider 



the first value of M(t)/t at the lowest value of t |I14I1 . We will 
use this latter estimate for deriving the drift and diffusion coef- 
ficients, underlying the evolution of NO2 concentration in Lis- 
bon. 

Within this framework, we consider the two-dimensional 
system of NO2 concentrations X = (xi , X2) describing the fluc- 
tuations at the stations of Chelas and Avenida da Liberdade, see 
Fig- [I] In order to comply with a Langevin process, as defined 
in Eq. we first verify that both data sets exhibit Markovian 
properties, which we show next for component xi only, for sake 
of clarity. For X2 the results are similar 

As Fig.[3]indicates, the conditional moments of the time se- 
ries show no evidence of measurement noise as t approaches 
zero Hell: M(iVt do not diverge when r ^ 0. This is true, 
both before and after de trending. For r smaller than a limiting 
value Tf, some oscillations are observed in the case without de- 
trending, although they have no impact on the estimate of the 
corresponding Kramers-Moyal coefficients, as compared to our 
method of using the value at t = 1 as estimate. For details see 
Ref. Hi . 

The resulting components of the drift and diffusion coeffi- 
cients are plotted in Fig. |4] As one sees all surfaces are ade- 
quately fitted by a quadratic polynomial 



P(Xl,X2) 



2 2 

■ aiX| -I- 02X2 + 03X1X2 -I- 04X1 + 05X2 + flg , 



(10) 



where p denotes the drift and diffusion components, D:'* and 



(2) 

D-j respectively, and the coefficients a, are computed from a 
least-square procedure on the drift and diffusion components as 
functions of the detrended variables xi and X2. 
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Figure 5: Contour plots of conditional probabilities (solid curves) and con- 
ditional two-point probabilities (dashed curves) computed from the detrended 
time series .vi and X2 with r2 = 1 hour for (a) rj = 2 hours and (b) 73 = 10 
hours. The con'esponding cuts through contour planes, indicated by the hori- 
zontal dashed Unes, are shown in (c) and (d) with a good matching between the 
respective one-point and two-point conditional probabilities. The distributions 
were computed with 13 bins for each variable using a sample of 10^ data points. 



4. Analysis of Markov properties for the series of NO2 con- 
centrations 

The Markovian nature of the variable xi can be investigated 
by considering the differences between the conditional one- 



point probability p{x"^\t\x'^' ,t - T2) and the conditional two- 
point probability p{x^^\ t\x^^\ t - T2; xf\ t - T3). If the process 
is Markovian on time scales larger than T2, then these prob- 
ability distributions should not differ significantlv lfioll for any 
choice of T3. Indeed, as can be seen from Fig.|5] the Markovian 
properties seem to be fulfilled for T2 = Ih and both T3 = 2h 
and T3 = lOh . We therefore observe strong indications that the 
process is Markovian already at the sampling rate of the data 
points of Ih and for time lags longer than Ih. 

Further, it is also necessary to check the Gaussian nature of 
the stochastic force F and ascertain it indeed obeys Eqs. (|4]i. Us- 
ing the measured time series and the estimated ICM-coefficients, 
the noise r(f) can be reconstructed from a numerical discretiza- 
tion of Eq. (O solved with respect to F[25], namely solving 



v(2) 



g-'(X)(X-h(X)), 



(11) 



where X = X(f -H 1) - X(f) and h(X) and g(X) are evaluated at 
X = X(f). 

The resulting noise is analyzed with respect to its autocorre- 
lation, shown in Fig. |6] the autocorrelation decays to zero for 
the very first values of t, which strongly supports to treat Fi and 
F2 as a white, 5-correlated noise source. 

For ascertaining the Gaussian nature of the stochastic sources 
we plot in the inset of Fig. |6]the PDF of the reconstructed noise 
time series Fi and F2 (solid lines) against a Gaussian distribu- 
tion (dashed lines). 



Figure 6: Autocorrelation of the reconstructed dynamical noises T\,Ti 
(stochastic fluctuation), indicating that they are r5-correlated. The inset shows 
the probability density function (PDF) of the reconstructed noise normalized to 
variance 1 (lines) and a normal distribution for comparison (dashed line). 



As one sees from the inset, in the range comprising over 
95% of the Gaussian noise, the distributions for the stochastic 
sources are well approximated by a Gaussian distribution. We 
find it reasonable to assume, therefore, that the data series can 
be approximated sufficiently well by a Fokker-Planck equation. 
The deviations observed for the extreme values, are common in 
the analysis of long-term field measurements, showing tails for 
close to exponential decay. 

From the tests described in this section one may satisfacto- 
rily take the series xi and X2 as a set described by two coupled 
Langevin Equations, Eq. (O with K -2. Next we derive these 
equations from the sets of measurements x\ and X2. 



5. Deriving optimal variables: eigensystem for NO2 mea- 
surements at different stations 

Having successfully determined the drift and diffusion 
constants describing the respective deterministic forcing and 
stochastic fluctuations of the system of NO2 concentration mea- 
surements, we now determine the eigensystem of the diffusion 
matrices and investigate its principal directions. This procedure 
was described in detail in 11911 and was previously applied to a 
two-dimensional sub-critical bifurcation t26] and to the analysis 
of human movement ll27ll . It will be briefly outlined here, for K 
variables. 

Diffusion matrices are numerically estimated on a mesh of 
points in phase space, as shown for example in Fig. |4j;-e. Then 
at each mesh point the K eigenvalues and corresponding eigen- 
vectors of the estimated matrices are calculated. The diffusion 
matrices contain information about the stochastic fluctuations 
acting on the system and we use the local eigensystems of the 
matrix for a further characterization of these forces. In par- 
ticular, a vanishing eigenvalue indicates that the corresponding 
stochastic force may be neglected. 
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We are looking for a transform of the original coordinates 
X = {xj} into new ones X = {i,}, such that the new coordinates 
are aligned in the directions of the eigenvectors of the diffusion 
matrix in each mesh point, i. e. the principal direction in which 
the diffusion matrix is diagonal. Diagonalizing the diffusion 
matrix decouples the stochastic contribution in the set of vari- 
ables, and if the eigenvalues in the transformed coordinates are 
significantly different, we are able to restrict our investigation 
to the coordinates with lower stochasticity. 

We therefore look for a two-times continuously differentiable 
function F with 

X = F(X,f), (12) 

for which||23, 24] the deterministic and stochastic parts in the 
Langevin systems of equations, transform respectively as lllQIl 



k=l 



dxk 



l/>(X) 



|]§*>(X) 



A=l 



dxk 



(13) 
(14) 



where the second equation reads g(X) - J(X)g(X), with J(X) 
the Jacobian of our transformation F. For reasons of clarity in 
the following we do not explicitly notate the dependence on X 
andX. 

The eigenvectors of matrices g with coordinates in 
local bases e,, can be incorporated in matrices U = 
[ui U2 ... vlk\- Defining U as U = J^U one then obtains 
(see Eq. (fT4l i) 



U^g^gU. 



(15) 



By definition the inverse transform F '(X) is chosen such 
that the normalized eigenvectors are given by 



1 d¥- 

Sk dxk 



with 



Sk = 



dF- 



dxk 



(16) 



(17) 



i.e. the respective square sum of the columns in the Jaco- 
bian of the inverse transform. Taking into account this scal- 
ing factor the eigenvalues in the new coordinate system can be 
calculatedin. 



D<2' = diag 



p2 



A2 
p2 



,.2 



(18) 



where Ai {i - \, . . . ,K) are the eigenvalues of the diagonalized 

f2) 

matrix D\ . 

diag 

In general, the eigenvalues of the diffusion matrix indicate 
the amplitude of the stochastic force and the corresponding 
eigenvector indicates the direction towards which such force 
acts. These directions can be regarded as principal axes of the 
underlying stochastic dvnamics lll9ll . 



In particular the vector field yielding at each mesh point the 
eigenvector associated to the smallest eigenvalue of matrix g 
defines the paths in phase space towards which the fluctuations 
are minimal. If this eigenvalue is very small compared to all 
the other, the corresponding stochastic force can be neglected 
and the system can be assumed to have only K - \ independent 
stochastic forces, reducing the number of stochastic variables 
in the system. 

Notice however that, whereas in a Cartesian coordinate sys- 
tem the eigenvalues are strictly related to the amplitude of dif- 
fusion in the corresponding eigenvectors, a nonlinear transfor- 
mation usually changes the metric|[l9, 28 1. In the transformed 
system the direction of the maximal eigenvalue is not neces- 
sarily the direction with the highest diffusion. This disparity 
is accounted for by the factor s, above. A much more simple 
case occurs when the eigenvectors are parallel (or almost) to a 
fixed direction, meaning that the eigendirections at each point 
in phase space are the same but rotated by a constant angle. 
Next we address such situation. 

6. Transform of NO2 concentrations to the stochastic 
eigendirections 

In this section we apply the procedure described previously 
to the two series of NO2 concentration, x\ and xi in Chelas and 
Avenida da Liberdade, shown in Fig.|7^. The joint PDF of both 
concentrations x\ is X2 is plotted in Fig.|7j5 showing the region 
in phase-space most visited by the bivariate series {x\,X2). A 
plot of the eigenvectors of D*^' in Fig.|7}; suggests that a contin- 
uous and smooth description of the corresponding sorted eigen- 
values exists. Here we place at each grid point of the phase 
space one ellipsoid whose major and minor axis are given by 
the (non-normalized) eigenvectors associated to the largest and 
smallest eigenvalue respectively. Since the two eigenvalues are 
different, the eigenvector corresponding to the lower eigenvalue 
describes the direction of minimum stochasticity. The eigen- 
vectors and eigenvalues of the diffusion matrix give locally the 
principal directions of stochastic fluctuations (diffusion). 

In general, from a plot as the one in Fig.|7}; it is possible to de- 
rive numerically the variable transformation in Eq. (fT2l i: at each 
grid point one determines the angle between the "largest" eigen- 
vector and the positive horizontal axis. Figure |7}1 shows the an- 
gle 4> for the bivariate series {x\,X2). Rotating each ellipsoid 
separately by the respective (^-angle aligns the largest eigenvec- 
tor along the horizontal direction and the smallest eigenvector 
along the vertical direction, yielding the two new (transformed) 
variables x\ and X2- This angle can be derived at each grid point 
from the corresponding diffusion g components namely 



tan 20 - 



gn -gii 



(19) 



The angle or its absolute value quantifies the relative off- 
diagonal contribution that describes the coupling of the noise 
terms by the diffusion matrix. 

In general, what does such a transformation add to our un- 
derstanding about the system? First, by definition the transfor- 
mation decouples independent stochastic forces in the system. 
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Figure 7: (a) The original time-series xi and X2 and (b) their joint probability 
density function. In (c) we plot the two (orthogonal) eigenvectors defining the 
semi-axis of an ellipsoid. The major semi-axis is associated with the largest 
eigenvalue and correspondingly the minor semi-axis with the smallest. For 
each grid point in phase space {xi,X2) we plot (d) the angle <l> between the 
eigenvector associated to the largest eigenvalue and the positive jr-semi-axis. 
As one sees from (c) and (d) the angle (fi does not show significant disparity in 
its values within the considered range, indicating a strong uncoupling between 
the two stations (see text). Thus, a global rotation of the axis may be considered 
(see Fig. [8). 



The original (detrended) pair of variables as well as the trans- 
formed pair of variable obey Eq. (O, with one important dif- 
ference: the transformed pair of variables are such that each 
variable has a stochastic contribution governed by one indepen- 
dent stochastic force alone. In other words g(xi,X2) is diag- 
onal. For the original pair of variables the stochastic contri- 
bution mixes both independent stochastic forces. Second, in 
a reference frame where the two independent stochastic forces 
are decoupled, their minimum and maximum magnitude reach 
the largest difference between them. In other words, one aligns 
the major and minor axis of the "diffusion ellipsoids" shown 
in Fig. I?};. In the particular case when one of the magnitudes 
is much smaller than other one, one of the variables can even 
be disregarded as a stochastic variable, reducing the number of 
stochastic variables describing the system. A generalization to 
K variables is straightforward. 

The value of {(p) shown in Tab.[T]is the (jci, X2)-averaged an- 
gle {(p) ^ 0.407r ~ 7t/2 indicated at the scale of|7}i. Further, one 
inspection of Fig.|7}; and|7}l enables the observation that in our 
present case the 0-angle remains approximately constant at any 
grid point. Similar observations were made for the other pairs 
of stations in Lisbon (not shown). Consequently, we may con- 
clude that for our set of stations a global rotation is enough to 
align the "diffusion ellipsoids". For the stations in Chelas and 
Avenida da Liberdade, Figure[8^ shows the result obtained after 
performing a global rotation by the median median(0) ^ 0.43;7r. 
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Table 1 : Characterizing dift'erent pairs of variables: the oiiginal pair of mea- 
sures y, the detrended pair of variables x, the transformed pair x, and, for com- 
parison, the pair x transformed according to the simpler rules in 1221 . For each 
vaiiable or pair of variables (/ = 1, 2) we show the mean (); and standard devia- 
tion cr, of their distribution of observed values together with the rotation angle 
averaged over phase space, as well as the average coefficients Q and Rj for 
evaluating their stochastic and deterministic contributions. See Eqs. U9t , )2Qt 
and i2l\ . The coiTelation coefficient between both variables is also given in 
each case (see text). 



In Fig. [SJ) both eigenvalues /I, are plotted, corresponding 
to the length of the major and minor axis of the diffusion el- 
lipsoids. While there is a significant difference between both 
eigenvalues, ~ 2.5 A„,i„ as shown in Fig.|8f, they are of the 
same order of magnitude. Such observation indicates the pres- 
ence of two independent stochastic forces driving the bivariate 
signal {xi,X2). 

The stochastic contribution for each variables of the pair 
(x\,X2) obeying Eq. ^ can be compared through one param- 
eter Q defined at each point x in phase-space as 



^ ~ 1 / \ , 2 / \ 



(20) 



where one orders the rows of matrix g to guarantee Q < I, 
i.e. variable xi is chosen as the one having lower stochastic 
contribution. When Q = I both stochastic contributions are 
equal. When 2 1 one stochastic contribution can be ne- 
glected, reducing by one the number of stochastic contributions 
in the system. For an arbitrary number of stochastic variables, 
the generalization of Eq. ( |20] | is strai ghtforward lfl9l . 

Table [T] shows the value of coefficient Q for the set of mea- 
surements y, for the detrended variables x and for the trans- 
formed detrended variables x. The coefficient is averaged over 
the sample of points in the corresponding phase space. For 
y and x the smallest stochastic contribution has a magnitude 
of approximately 70% of the largest one, while for the trans- 
formed variables it decreases more than 2%. This magnitude is 
not small enough to permit neglecting one variable. We con- 
sider this finding the central result of this letter; before trans- 
formation the pair of detrended variables include already two 
independent stochastic forces of the same order of magnitude. 

One note is however important to stress at this point. The 
method applied here to empirical data deals with a transforma- 
tion that operates on the diffusion matrix alone. No constraints 
related to the drift functions, hi and I12 are considered. To eval- 
uate the predictability of each variable / one needs to compare 
the total amplitude of the stochastic term with the deterministic 
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Figure 8: Transformed variables thi'ough a global rotation of the xi and X2 axis (check Fig. |7) into new variables, .vi and X2. (a) Eigenvectors of the transformed 
variables and (b) its eigenvalues derived for the diffusion matrix of the transformed variables, together with quantities to evaluate some underlying properties of the 
system, namely (c) the rotation angle (see Eq. (19)), (d) the asymmetry of the stochastic influence at each variable, given by Q in Eq. j20K (e) the deterministic 
coefficient R in Eq. )2U and (f) the quotient between the largest /l„,„v and the smallest /i„„„ for each grid point in the transformed phase position. For each property 
the average value is showed with horizontal surfaces and an explicit indication at the vertical axis. 



term, namely 



7. Comparison with standard methodologies 



Rj(x) = 



(21) 



Such expression is also straightforwardly extended to K vari- 
ables. The larger/?, the more predictable the variable ; may be, 
i.e. the smaller the stochastic overall contribution is compared 
to the deterministic part governing the evolution of the variable. 
In our present case, as given in Tab. \T\ while the detrending 
y ^ X of our measurements increases the predictability of the 
non-periodic modes in time, the global rotation has no major ef- 
fect: both coefficients R, maintain the same order of magnitude 
after transform. 

The correlation coefficient /i between both stations is also 
given in Tab. [1] While detrending has no significant effect on 
the correlation, the transform x, — » i, indeed decreases its ab- 
solute value. 

Figure [8};, [8]i and [8^ illustrate the numerical result of each 
property, 0, Q and R for the transformed variables. Similar to 
such variables is the quotient between the maximum and mini- 
mum eigenvalues, shown in Fig.[8f. Similar plots are obtained 
for the other possible pairs of stations. 



In this Section we first address the question of how good the 
coordinate transform derived above is compared to other, pos- 
sibly simpler transforms. 

For example, we may consider a transform to coordinates 
which describe the mean value and difference between the two 
measured time series, e.g. 



Xi + X2 



xi 



Xi 



(22) 



This choice is the simplest one for two variables, one describ- 
ing the total amount x\ + X2 and another describing the relative 
amount xi - x^- For such choice of variables we obtain a value 
of 2 = 0.59, which is essentially the same as for our "opti- 
mized" variables (see Tab. [TJ, The absolute value of the an- 
gle (101) is however considerably larger than for our optimized 
transform, as is the correlation coefficient between the time se- 
ries, meaning that this simple transform fails to decouple the 
noise sources. The drift-diffusion quotients yield /?i - 0.17 and 
Ri - 0.24, showing again no better predictability in comparison 
with the original variables. 

In our case we saw that the eigendirections do not depend 
much on the detrended variables x\ and xi, which implies that 
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they are functionally decoupled. However, sometimes it is nec- 
essary to consider a proper scaling of the variables lll9ll . In such 
cases, we find it advisable to use a more general transform to 
generalized polar coordinates given by 



^ctan(|^) + e 



(23) 



where in general the radial and angle variables, and 9g, are 
functions of the detrended variables x\ and X2. This approach 
has the advantage that the inverse transform F"' is given by the 
simple form 





f 







i[rjCos(^?, 



(24) 



In addition, the metric factors introduced by the polar transform 
can result in a more pronounced separation of the eigenvalues 
in the transformed coordinates. 

The other question addressed within this section is the com- 
parison between methods applied to choose the most appro- 
priate inputs. Theoretically, any set of input data can be fed 
into a model for training and evaluation. However, the num- 
ber of possible variables to be used and the number of ways 
they can be presented is too diverse to test all possible com- 
binations. A number of statistical methods can be appUed in 
order to choose the most appropriate set of predictors or inputs. 
Examples are, among other, stepwise regression, PCA, cluster 
analysis and ARIMA. For details, see Ref. ll22ll and references 
therein. Such pre-processing procedures, reduce the number 
of input variables into the models, thus eliminating the redun- 
dant information. In these standard procedures, the selection 
of variables is usually made independently for each monitoring 
station. 

Another possible way to tackle redundancy is the pre- 
processing of data consisting of the computation of backward 
stepwise regressions (BSR) conducted between a target variable 
and all the other data sets. Based on the available common pe- 
riod data sets, one constructs a collection of records, composing 
the input vector, which includes the meteorological variables, 
air pollutant concentrations, etc, and together with it assumes 
the corresponding target, which in our case is the atmospheric 
concentration of a certain pollutant. Subsequently, one retains 
the smallest subset of statistically significant variables to pre- 
dict a certain pollutant concentration automatically at a given 
monitoring station. In addition, BSR allows the determination 
of the best time lags for each input variable, typically daily and 
weekly cycles. 

The referred techniques also allow the comparison between 
the original data sets and surrogate data sets including only the 
stochastic component. The stochastic component may be deter- 
mined through a rough approximation of a mathematical func- 
tion (e.g., sin x), or, for example, by the presented framework. 
After the selection of variables and the determination of cyclic 



and stochastic behaviors on each time series, linear and non- 
linear models can be applied in order to model air pollution in 
each monitoring station. The forecasting capabilities of the dif- 
ferent approaches can then be compared. Such models are also 
applied to each decoupled time-series in order to predict next 
days air quality at each monitoring station. The applications 
of this framework, however, allows to determine the stochastic 
component on a efficient manner, enhancing air quality predic- 
tions. 

8. Discussion and Conclusion 

In this paper, we investigated the stochastic properties of a 
set of two simultaneous series, obtained by introducing a proper 
detrending of NO2 measurements, which is able to remove pe- 
riodic modes in the series. We focused in the measurements at 
two different stations out from a set of 22 stations in Lisbon. 

Based on validity tests we assumed, that the time series af- 
ter detrending were properly modeled by a system of Langevin 
equations. The validity of this assumption is discussed in sec- 
tion[3] showing that the data sets obey the Markov property to a 
sufllcient extent. The stochastic fluctuations show good resem- 
blance with 5-correlated Gaussian noise. 

Calculating the eigenvalues of the diff'usion matrices, we 
found a transform that leads to a description in which the dif- 
fusion matrices are diagonal. Since the transformed variables 
are derived directly from the transformation that diagonalizes 
the diffusion matrices, they correspond to the orthogonal direc- 
tions in phase space in which fluctuations are stronger (larger 
eigenvalue) and weaker (smaller eigenvalue) respectively. 

Comparison between original and transformed variables 
showed that the two detrended variables are driven by stochas- 
tic forces almost decoupled from each other, showing an almost 
constant rotation angle of the "diff^usion ellipsoid" at each point 
of phase space. Further, both stochastic sources have ampli- 
tudes of the order of the deterministic terms, indicating a short 
horizon of predictability. This procedure worked out well for 
the NO2 data, since the transformation of variables resulted in 
decoupling the diff'usion components in the new coordinates. 
Other transformation could be considered. For instance, we dis- 
cussed how this approach could be applied for other data sets 
in which the diffusion ellipsoids do not align in phase space, 
but instead depend in non trivial functional of the variables. In 
this case, the transform maps the detrended variables into two 
polar-like coordinates. 

One question that should be addressed in a forthcoming study 
is to present a systematic overview on all pairs of stations stud- 
ied by us in this scope but not shown thoroughly, since it was 
out of our main purposes. Doing that one would be able to com- 
pare in detail the results obtained through the method applied 
in this paper with standard methods used for forecasting NO2 
concentration at a specific spot in the city of Lisbon. 
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